
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# DISCLAIMER AND GENERAL INFORMATION
#
# File: SI5_robustness.R
# Purpose: Produces results presented in section SI5
# Date: June 2025
# Data: pulled through 00_data_prep.R
#
# See 00_data_prep.R for technical disclaimer on R versions
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Load data
source("00_data_prep.R")

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (A) Robustness: No control variables (Figure A5) ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Estimate model
m2 <- multinom(factor(effect_cc_1st) ~ sample, data=full)

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p5)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p5)==FALSE & p5=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI5_robust_wocontrols.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (B) Robustness: Social class controls (Figure A6) ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Define control variables to proxy social class
ctrls2 <- paste(c("factor(educ_cat)", "factor(income)", "factor(empl_sector)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls2, sep="+"), data=full)

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p5)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p5)==FALSE & p5=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI5_robust_social_class.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (C) Robustness: Exclude general population sample (Figure A7) ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full[full$sample!="General population",])

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,3,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 2,
                         grepl("Wales",hypothesis) ~ 3),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt, aes(x=estimate, y=sample, group=group, 
                         color=p5)) +
  # geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
  #           fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8), key_glyph = "pointrange") +
  geom_text(aes(label = ifelse(is.na(p5)==FALSE & p5=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI5_robust_drop_genpop.pdf", width=6, height=4)




# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (D) Robustness: New reference category (Figure A8) ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

table(full$residence, full$survey, useNA="ifany")
table(full$gor_england[full$survey=="Prolific"])

# Create new baseline category in general population sample
# Exclude respondents in general population data consisting of England minus North East, North West, and Yorkshire GORs
smb <- full[!(full$survey=="Prolific" & (full$gor_england=="North East" | 
                                           full$gor_england=="North West" | 
                                           full$gor_england=="Yorkshire and Humber")),]


# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=smb)

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))


# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
  filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
           hypothesis=="(1 Wales) - (1 Northern England)" |   
           hypothesis=="(2 Scotland) - (2 Northern England)" |
           hypothesis=="(2 Wales) - (2 Northern England)" |
           hypothesis=="(3 Scotland) - (3 Northern England)" |
           hypothesis=="(3 Wales) - (3 Northern England)" |
           hypothesis=="(4 Scotland) - (4 Northern England)" |
           hypothesis=="(4 Wales) - (4 Northern England)") %>%
  mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
                         grepl("Wales",hypothesis) ~ 4),
         group=case_when(grepl("1",hypothesis) ~ 1,
                         grepl("2",hypothesis) ~ 2,
                         grepl("3",hypothesis) ~ 3,
                         grepl("4",hypothesis) ~ 4),
         p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
         p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
  select(rowid,group,hypothesis,p5,p10) %>%
  right_join(d2,by = join_by("rowid","group")) %>%
  arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                color=p5)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p5)==FALSE & p5=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/SI5_robust_ref_category.pdf", width=6, height=4)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#                         END OF FILE
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

